function [gra_x,gra_y] = grad(I)
[m n] = size(I);
gra_x = zeros(m,n);
gra_y = zeros(n,n);
gra_x(1,1) = -2*I(2,1) - I(2,2);
gra_x(1,n) = -I(2,n-1) - 2*I(2,n);
gra_x(m,1) = 2*I(m-1,1) + I(m-1,2);
gra_x(m,m) = I(m-1,n-1) + 2*I(m-1,n);
gra_y(1,1) = -2*I(1,2) - I(2,2);
gra_y(1,n) = 2*I(1,n-1)+I(2,n-1);
gra_y(m,1) = -I(m-1,2)-2*I(m,2);
gra_y(m,n) = I(m-1,n-1)+2*I(m,n-1);
for j=2:n-1
    gra_x(1,j) = -I(2,j-1)-2*I(2,j)-I(2,j+1);
    gra_x(m,j) = I(m-1,j-1)+2*I(m-1,j) + I(m-1,j+1);
    gra_y(1,j) = 2*I(1,j-1)+I(2,j-1)-2*I(1,j+1)-I(2,j+1);
    gra_y(m,j) = I(m-1,j-1)+2*I(m,j-1)-2*I(m-1,j+1)-I(m,j+1);
end
for i=2:m-1
    gra_x(i,1) = 2*I(i-1,1)+I(i-1,2)-2*I(i+1,1)-I(i+1,2);
    gra_x(i,n) = I(i-1,n-1)+2*I(i-1,n)-I(i+1,n-1)-2*I(i+1,n);
    gra_y(i,1) = -I(i-1,2)-2*I(i,2)-I(i+1,2);
    gra_y(i,n) = I(i-1,n-1)+2*I(i,n-1)+I(i+1,n-1);
end
for i=2:m-1
    for j=2:n-1
        gra_x(i,j) = I(i-1,j-1)+2*I(i-1,j)+I(i-1,j+1)-I(i+1,j-1)-2*I(i+1,j)-I(i+1,j+1);
        gra_y(i,j) = I(i-1,j-1)+2*I(i,j-1)+I(i+1,j-1)-I(i-1,j+1)-2*I(i,j+1)-I(i+1,j+1);
    end
end